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The equations describing nonradial adiabatic oscillations of differentially rotating relativistic stars 
are derived in relativistic slow rotation approximation. The differentially rotating configuration is 
described by a perturbative version of the relativistic j-constant rotation law. Focusing on the 
oscillation properties of the stellar fluid, the adiabatic nonradial perturbations are studied in the 
Cowling approximation with a system of five partial differential equations. In these equations, 
differential rotation introduces new coupling terms between the perturbative quantites with respect 
to the uniformly rotating stars. In particular, we investigate the axisymmetric and barotropic 
oscillations and compare their spectral properties with those obtained in nonlinear hydrodynamical 
studies. The perturbative description of the differentially rotating background and the oscillation 
spectrum agree within a few percent with those of the nonlinear studies. 
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I. INTRODUCTION 



Differential rotation can appear in many phases of the stellar evolution of a neutron star, such as in protoneutron 
stars 

0,1 

, in the massive remnant of binary neutron star mergers [H, H, or as a results of stellar oscillations (r-modes) 
that may drive the star into differential rotation via nonliner effects [a, [(| 0, H| ■ The differentially rotating phase can 
last at most some seconds depending on the dissipative mechanism that drive the star to uniform rotation, such as 
magnetic bracking [t| [l(| and turbulent motion [ll|, [l2j • 

However, this very short period can appear during the most violent phases of neutron star's life, as in the case of a 
core collapse or binary merging. This is exactly the time when we expect to get the strongest emission in gravitational 
waves. Since the ground based detectors reached sensitivities which allow the detection of gravitational wave signals 
from oscillating or unstable neutron stars the exact frequencies of the emitted waves are urgently needed. The study 
of linear pulsations and stability of differentially rotating stellar models started more than three decades ago both in 
Newtonian and relativistic gravity (l3l. Il4l. [Tol [ltl . The spectrum of differentially rotating stars in Newtonian theory 
has been studied in the frequency domain by [17| . 

Recently, the effects of differential rotation on the dynamical and secular instabilities of the r and f-modes gained 
more attention. In Newtonian theory, the r-mode spectrum has been studied in [18l Il9l . |20( , while the f-mode and the 
secular stability limits have been investigated in [21| . In stars that rotate with a high degree of differential rotation 
(51 c /51 s « 10), an m = 2 dynamical instability can appear even for considerably low rotation rates (T/|VF| ~ O(10 -2 )) 
as suggested in [22|, [23[ • Note that T is the rotational kinetic energy, W the gravitational potential energy, Q c and 
Q s the angular velocity at the center and at the equatorial surface, respectively. In addition, an m = 1 dynamical 
instability has been identified for high degrees of differential rotation and soft equations of state 0, Hf|- More 
recently, it has been reported the discovery of to = 1 and m = 2 dynamical instability even for stiff equations of 
state [26]. Studies based on linear analysis [U [28| suggest that low T/|W| instabilities might be triggered when 
the corotation points of the unstable modes fall within the differentially rotating structure of the star. It is then 
evident that differentially rotating compact stars got a renewed interest and there are several recent fully general 
relativistic studies related to the issue. For example, slowly and differentially rotating magnetized neutron stars have 
been studied in [29, 30j, the collapse and black hole formation of magnetized and differentially rotating neutron stars 
in [3lj and the frequencies for nonlinear axisymmetric pulsations of differentially rotating stars in the Cowling and 
conformal flatness approximations in [32l . |33| . 

As a further contribution to the direction of understanding stellar differential rotation we present here the equations 
describing a slowly and differentially rotating relativistic star. In our study we use a perturbative analysis and keep 
terms up to first order with respect to stellar rotation. In this first perturbative approach, we neglect all the spacetime 
perturbations i.e. we use the so called Cowling approximation since we want to focus on the behavior of the stellar fluid. 
The slow rotation approximation has been extensively used in the study of stellar perturbations both in Newtonian 
and in general relativistic approach [34l l35l . l36l |37| . Here we study the spectrum of axisymmetric perturbations in 
order to test the accuracy of our approximation technique against published results by nonlinear numerical codes 
[HI, l33l ]. The next steps will be to estimate the oscillation spectrum of the non-axisymmetric perturbations and 
furthermore to proceed in studying the low T/|W| dynamical instability. 

The structure of the paper is as follows. In section [II] we construct the differentially rotating stellar models, and 
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test their accuracy against the equilibrium configurations obtained with non-linear codes. In section III, we derive 
the perturbation equations in the general case of non-barotropic and non-axisymmetric perturbations, which are 
numerically solved in section llVl for the axisymmetric and barotropic case. In the section [V] conclusions are drawn 
and the future applications are proposed. The appendix is organized in four sections. Section [X] is dedicated to 
the analytical expressions of the harmonic components of the j-constant rotation law, while all the coefficients of 
the perturbed conservation equations are given in section [Bj The operators that couple perturbations with different 
harmonic indices are introduced in section [Cj Finally, in the last section [D] we discuss the j-constant rotation law in 
isotropic coordinates. 

Throughout the paper we use geometrical units c = G = 1. We use a prime ( ' ) to denote derivatives with respect 
to the radial coordinate r and the overdot ( ' ) for derivatives with respect to the time coordinate t. 



II. EQUILIBRIUM CONFIGURATION 



The axially symmetric spacetime of a slowly and differentially rotating star can be described by the following line 
element: 



ds 2 



^dt 2 + e 2X dr 2 



2 tor 2 sin 2 6 dt < 



(1) 



where A are scalar fields that depend on the radial coordinate r. The metric function w, which describes the 
dragging of the reference frames due to the stellar rotation, is a function of the both r and 8 coordinates. In this 
paper, we assume that the stellar interior is described by the perfect fluid energy-momentum tensor: 



T al3 = (e+p) u a u f3 + pg af) , 
where e and p denote the total energy density and the pressure respectively, and u a is the fluid velocity 

u a = ( e -", 0,0, fie"") , 



(2) 



(3) 



where = fi (r, 9) is the angular velocity of the star as measured by an observer at infinity. 

By providing an equation of state (EoS) p = p (e), a rotating equilibrium configuration can be determined by solving 
the Tolman-Oppenheimer-Volkoff (TOV) equations together with the following equation for the metric function uj: 



Ait (e + p) re 



2 A 



167r (e + p) 



1(1 + 1) 



„2A, 



-16tt (e + p) e 2X n 



(4) 



A. Rotation Law 



In general the differential profile of the angular velocity Q in not known. Still approximate analytic differential 
rotating laws have been adopted and used both for Newtonian and General Relativistic configurations (see [38| 
and reference therein). Among these laws, the j-constant rotation law [3^,|43| guarantees Rayleigh's stability against 
axisymmetric perturbations for inviscid fluids, and seems to describe quite well not only typically differentially rotating 
stars but even the remnant rotational configurations that arise from the mergin g of hypermassive binary systems [H . 

The relativistic j-constant rotation law can be perturbatively expanded as in [4l| : 

A 1 + e lv r A sin 

where Sl c denotes the angular velocity at the rotation axis, while A is the parameter that describes the degree of 
differential rotation of the star. For high values of A (e.g. A ~ 500 km), the j-constant rotation law ([5]) tends to a 
uniformly rotating configuration, where £1 — > fl c . 

In the slow rotation approximation, rotation is treated as a perturbation of a spherically symmetric background 
spacetime. Therefore, the tensor harmonic expansion enables us to separate the angular dependence of any pertur- 
bative fields from the temporal and radial parts. In this way the metric function lu and the angular velocity f2 of the 
fluid can then be written as follows: 

u(r,9) = -£^)-^§, (6) 



sin6> 89 ' 
1 BP 

Q(r,6) = -^tir)— y, (7) 
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where P/(cos0) are the Legendre's polynomials. The Hi and u>i components in equations ©-([7]) are then determined 
via the orthogonality relations of the form, 

where / = f(r,9) is a scalar field. In equation ([8]), dft denotes the volume element of the 2-sphere S 2 . It is worth 
mentioning that for a uniformly rotating star, the only non vanishing components of equations ([6|)-([7]) are the ones 
for I — 1 [42]. In this case, u> depend only on the radial coordinate r and O is of course constant. For differentially 
rotating stars, which are described by the j-constant rotation law, the components Qi can be determined by direct 
substitution of equation (O in ([8]). Then the resulting expression can be split in two parts: 

fi,=2jf +2?, (9) 

where Xf is the "nearly Newtonian" part and Xf- is the "relativistic" correction. The definition of these two terms is 
given by the following expressions: 

rN _ 21 + 1 Q c f A 2 sm9 dPi 



X N — c I dVu (10) 

1 l(l + l)4n J S 2 A 2 + e~ 2 »r 2 sin 2 9 86 ' y ' 



Z R^ + 1 e- 2 V f - sin 3 ^M) dPj 
1(1 + 1) 4tt J S 2 A 2 + e~ 2v r 2 sin 2 9 89 ' 



The parity of the j-constant rotation law implies that the non- vanishing components in (|10[) - (|11[) are the ones with odd 
I. The integral Xf for I = 1 and I = 3 has been already calculated in [4l|, and is analytically written in appendix fAI 
Moreover, the integration of the "relativistic" term Xf- requires an additional expansion of the metric variable u>. 
Therefore, by inserting equation © into (fTTj) one gets 

i' 

where are given by the following integral relation 

Jl -J s2 A 2 + e- 2 »r 2 sin 2 9 89 89' [l6 > 

The detailed forms of Jl , and J% are shown in appendix [XJ Our choice to retain only the I — 1 and / = 3 
terms for the expression (|13p will be justified later in section HVl were we show that the frequencies of the nonradial 
oscillations that we calculate using this approach are in very good agreement with the ones derived using the full 
nonlinear evolution of the fluid background [321 . |33| . 



B. Frame dragging in differential rotation 

The various components lui of the series expansion © which provides ui will be determined from the solution of 
equation ^ for the frame dragging. In practice, one substitutes equations © and ([7]) into equation and collects 
the terms of the series with the same index I. However, the term Xf- describing the relativistic correction, equation 
(|I2[) . introduces a coupling between the different harmonics, which is due to the presence of metric function uj in 
equation ([5]). As a result the terms uji will be determined as solutions of a system of coupled ordinary differential 
equations. In order to simplify the presentation of this system of differential equations we will write it schematically 
as: 

C(u)=S(Sl), (14) 

where £(■) represents the linear differential operator of the left hand side of equation (UJ), while S(-) denotes the linear 
operator of the source. When equation (fT4)) is expanded into vector harmonics, it reduces to the following set of 
equations 

'W-'W + 1PT?^M')' (15 > 

which applies to any value of I. The first two non- vanishing terms lo\ and uj^ obey the following coupled system of 
equations: 
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Model T c (ms) Q c (xlCT 2 ) km' 1 ft e (x!0~ 2 ) km" 1 

Bl 1.719 1.218 0.435 

B2 1.204 1.740 0.621 

B3 0.970 2.160 0.771 

B6 0.657 3.189 1.139 

B9 0.496 4.219 1.507 



TABLE I: Properties of the background models used in the paper. The quantities T c and Q c are respectively the period and 
the angular velocity at the rotational axis, while Q e represents the angular velocity at the equator. All stellar models have 
mass M = 1.40 M Q and radius R = 14.151 km. 



J. r 2 e -2^ (wiJi i) = S(l^)+^e- iv S(u a Ji) , (16) 
CM-^e-^S(^Ji) = Sil^ + ^e-^Si^Ji) , (17) 

which together with a set of appropriate boundary conditions at the stellar center and at infinity form a boundary 
value problem. Actually, the demand for regularity of the solutions at the stellar center and their decaying asymptotic 
behavior suggest the following approximate relations [4^ |: 

uji ~ r 1 - 1 , for r^0, (18) 
uji — r~ l ~ 2 , for r^oo. (19) 

The numerical method used for the solution of the above system of equations is described in the next paragraphs. 



C. Stellar Models 



In order to test the accuracy of the perturbative treatment of stellar differential rotation, which has been introduced 
in this paper, we will compare our results with those derived from the numerical solution of the nonlinear Einstein 
equations. Thus we will study a set of stellar models, the so-called B models, which have been already adopted in 
[32 . H|[ . These models describe barotropic and differentially rotating stars, where the rotating pattern is described 
by the j-constant rotation law. For simplicity, a relativistic barotropic EoS of the form 

p = Kp r , (20) 

e = P+fTT' ( 21 ) 

has been adopted, where p is the rest mass density and K and T the polytropic parameters. 

For a given EoS, an equilibrium configuration depends on three parameters: the central density p c , the stellar spin 
at the rotation axis fi c and the parameter A describing the degree of differential rotation. For the B models the 
polytropic parameters get the values K = 217.858 km 2 and T = 2 while the central rest mass density is fixed to the 
value p c = 5.87 x 10~ 4 km~ 2 . The solution of the TOV equations provide a stellar model with mass M — 1.40 Mq and 
radius R = 14.151 km. The angular velocity at the rotation axis, f2 c , for some of the B1-B9 models is listed in tablefl] 
In order to compare our results with those of nonlinear studies special care has been taken for the interpretation 
of the parameter A, due to the different coordinate systems that are commonly used in non-linear studies and our 
perturbative approximation. This issue is explained in detail in appendix [D] where we show that A should be set to 
the value of the isotropic stellar radius. 

As we mentioned earlier the two components u)\ and u>z of the metric function u will be the solutions of the 
boundary value problem described by the two differential equations p6[) - (|17[) together with the boundary conditions 
(fT5]) - (fT!)|) . By discretizing equations (fH))) - (fl~T)) using a second order scheme one can construct a tridiagonal linear 
system. This system can be solved using a standard LU decomposition where the boundary conditions (fT8)) and (fT9]) 
are appropriately implemented. The presence of coupling terms on the right hand side of the equations calls for an 
iterative treatment. As a first step a trial solution ujj 1 ' of (jTSJ) - ((TTJ) can be determined by neglecting the coupling 
source terms, i.e. = j\ = 0. This trial solution is then substituted only into the source terms and the system of 
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r[km] 



FIG. 1: The harmonic components fh (upper panel) and H3 (lower panel) of the angular velocity Q of the Bl stellar model 
are plotted, where O c = 1.221CP 2 km -1 and A — 12 km. The dashed lines display the contribution of the "nearly Newtonian 
term" given by equation (|10p . 




r[km] 

FIG. 2: The harmonic components ui (upper panel) and u>3 (lower panel) of the angular velocity u> of the Bl stellar model are 
plotted. The dashed lines denote the trial solutions of equation (|16[) . which have been determined without the coupling terms 
and J\ . The solid lines are instead the final solutions. 



equations (|16[) - (|17[) is solved. This procedure is repeated until the solution converges. In practice, only three or four 
iterations are needed. 

In figures[T]|2J we show the behavior of the angular velocity f2 and the metric function ui for the slowly rotating model 
Bl, see table HI In figure [1] it is noticeable the small contribution of the "relativistic corrections" X l R in the estimation 
the star's angular velocity. The two components of the metric function i.e. uj\ and UJ3 are plotted in figure [H together 
with initial trial solutions w ; tr . It is worth noticing an imperceptible correction to the initial trial solution for 1 = 1, 
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FIG. 3: The dependence of the two components of the angular velocity Q on the differential parameter A is shown in this figure 
for the stellar model Bl. 



and the small value of 0J3 with respect to u>i, which suggests that terms of the perturbative series I©-® with higher 
I will not affect the results. The dependence of f^i and SI3 on the the differential rotation parameter A is illustrated 
in figure [3] There one can notice that for increasing values of A the £l\ reaches its uniform rotation value f2 c , while 
1^3 decreases proportionally as expected. 

In order to estimate the accuracy of our approximate solutions for f2 and u), we have done a comparison with the 
profiles derived by the nonlinear numerical code RNS (HJ . The results from RNS and the perturbative solutions for 
the Bl-model are compared in figure|3]on the equatorial plane. The maximum deviation between the two estimations 
is around r — 6 km where the relative error is about 5.5 % . The accuracy could be further improved by taking into 
account the / = 5 contributions. However, the nonradial perturbation equations would become considerably more 
complicated due to the presence of extra coupling terms. As we show in the next section the approximation used 
here allows the estimation of the nonradial oscillations frequencies with an error smaller than 10 percent even for very 
rapidly rotating models, see table HTT1 



In this section, we derive the perturbation equations which describe the nonradial oscillations of a slowly and 
differentially rotating star. The rotation, as described in the previous section, is treated perturbatively up to the first 
order in f2 [45j. In this work we focus on the spectral properties of the fluid modes, therefore, we neglect all metric 
perturbations, and we use the so-called "Cowling approximation". For slowly rotating stars, this method enables us 
to determine the frequencies of nonradial pulsations with an accuracy to better than 15 % with respect to the non 
approximated approach. In addition, the accuracy improves for modes with higher harmonic index I and mode order 

We will consider adiabatic oscillations, where the pressure and density perturbations are related by the following 
expression: 



where £ r is the radial component of the fluid displacement vector £^ and T and Ti represent the adiabatic index of 
the background and perturbed configurations respectively. The adiabatic index T is defined by the following relation: 



III. PERTURBATION EQUATIONS 




(22) 



r = 



p + e dp 



(23) 



e de ' 



7 




2 4 6 8 10 12 14 

r[km] 



0.004 - 



0.003 - 
3 0.002 - 




FIG. 4: The profiles of f2 and u) derived from the perturbative treatment and the nonlinear RNS code compared in this figure. 
The solid lines represent the perturbative solutions and dashed lines the RNS quantities. In the upper panel the estimations 
of the angular velocity S7 are plotted, while in the lower panel the variable w is compared to the shift component — /3 , which 
have been determined with the RNS code. 



and sound speed by 



2 _ £i 
Cs r de ' 

In Cowling approximation, the displacement vector £ M is related to the 4-velocity perturbation by 



5u„ = 



The r component of this equation provides an evolution equation for £ r : 

(d t + nd^C = e- 2X+v 5u r . 



(24) 



(25) 



(26) 



The three scalar quantities describing the perturbations of the total energy density, pressure and the r component 
of the Lagrangian displacement can be expanded in spherical harmonics, 



5e = SeimYim , 
Sp = Spi m Yi m , 



e 



v 1 

r 



Amy 
s 1 " 



(27) 
(28) 

(29) 



The perturbation equations can be considerably simplified if one uses a new function H, which is related to the energy 
density and pressure perturbations by the following relations 



Spi m = Hi m (p + e) 
(p + e) 2 



de 



(Hl m — £,lm) 



(30) 
(31) 



For barotropic fluids, H becomes the perturbation of the relativistic enthalpy [47J • Velocity perturbations are also 
expanded into vector harmonic expansion and have the form: 
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Su r = -e v }^ui t i m Yi m , (32) 



! 111 



Sue = e 2. (j" 2 < im ~W ~ U3 ' lm ~80~ ^ffifl ) ' (33) 

Ira v / 

Su^ = -e»J2 (u2j m ^+u 3j7n sm6^) , (34) 



i in 



5u t = ~^l(r,9)Su^ = -e" V" V" (^.Zm^rr 1 + M 3. Zm sin6>^^- ] £W^ r ) ^-7; "l^r j ( 35 ) 



Z'm' Zr. 



notice that in the last expression equation has been used. 

Due to the parity of the tensor harmonics, the nonradial perturbations can be divided in two classes, the axial (odd- 
parity) and polar (even-parity) perturbations. The polar sector can be described by four perturbation functions: the 
perturbation of the relativistic enthalpy H(r, t), the two velocity perturbations u± (t, r), U2 (f, r), and the r component 
of the displacement vector, £ r (r, t). On the other hand, in the Cowling approximation, axial perturbations will be 
described only by the velocity perturbation function 1*3 (t, r) . Using the perturbed energy- momentum conservation 

equations 5 (T a p'^j = one can derive a system of coupled differential equations for Hi rn and the Uj,2 m which together 

with equation f|26[) form a complete linear system of evolution equations for the five unknown perturbation functions. 
The four components of the perturbed energy-momentum conservation equations get the form: 



[t] 


a (*) y. + 


R (i) 




8Y lm 
86 


[>■] 




~B {r) - 




dY lm 
86 



[0] 



0. 



9 (6) 9 

aij m + a 2 ,i m cos 9 + a 3 j m sin 9 + a\J m cos 9 sin I 

9 (G) 9 

bi,i m + h.im cos 9 + b 3t i m sin + b\ [ m cos 9 sin 9 



dY lm 
89 

1 8Y ln 



[<j>] bij m + b 2 ,im cos 9 + b 3t i m sin 2 9 + bf} m cos 9 sin 2 ■ 



01,2m + a 2,im cos + a 3 ,; m sin 2 9 + af] m cos 9 sin 2 I 



sin 
0*z 



. 8Yl m 



0. 



89 
1 Wh, 



41 + 41 Bin 2 < 



(36) 
(37) 

(38) 

sm9Y lm =0, (39) 



where all the coefficients Ai m , Bi m , Ci m , di t i m , bi t i rn and Cj,; m are analytically written in appendix IB1 Notice that in 
equations (|26p and (|36[) - (|39 [) . the background variables lu and Q are expanded up to the harmonic index I = 3. The 
angular dependence of the above system of equations is removed by performing angular integrations. As a result, we 
obtain the following system of coupled evolution equations for Hi m , ui,2 m U2,Zm, U3,im and £[ m : 
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Zmax F (kHz) Hi (kHz) 

2.687 4.551 

1 2.710 4.571 

2 2.712 4.575 

3 2.712 4.575 



TABLE II: Oscillation frequencies of the quasi-radial F and Hi modes for the Bl stelar model. With Z max we denote the 
maximum harmonic index that has been used in the estimation. It is apparent that there is a fast convergence towards the 
exact value for relative small values of Z max - 



A 



I'm 



4m + im B\ 



Aoi, 



D 



(t) 



1,1m 



Cf 2 B. 



2.1m 



>(r) 



32, In 



£-3 1 (l2,lr 



, r±2 R (r) 
+ L 1 B 2.lr, 



= 0, 



*-l ( - / l,im +Z -l °2 



y 2.lr. 



£ ±3 a W 

^2 M 4,Zm 



0, 



-3 u 4.Z?r 



+ C W 



2^1 (0) 



= 0, 



±2 



-2,ir. 



(40) 
(41) 

(42) 



A6i„ 



0.2,1m 



+ Cf 1 b2Ar. 



- 2Cf 1 a,3j r , 



±2 (fl) 



^2 w 4,im 



±2„M 



2 r ±l,(S) 



m z £f L b 



4,lm 



C l,Zm 



/•±3 (40 n 
*-3 c 2,Zm U - 



(43) 



& m -I- im (fiio + 6S1 30 ) &r, 



1 Wj v'u\.l m + — O. 30 £f 2 ^l m 



(44) 



where A = 1(1 + 1). The operators £f J couple perturbations with different harmonic indeces I and they are defined 
in appendix [C] The final equations can be written in a more compact form as, 

Pirn + im [P lm + Ai±i^ rn + Pl±2, m J + M±.l,m + ^Z±3,m = , (45) 
A irn + Vn \ Al m + Pl±l, m + Ai ± 2,m^j + Pl±l,m + Pl±3,m = , (46) 

where Pi m and Ai m represent polar and axial perturbation functions respectively. The tilded quantities denote the 
extra coupling terms introduced by the differentially rotating background. It is worth noticing that apart from the 
I ± 1 couplings existing in the case of uniform rotation we have extra couplings I ± 2 and I ± 3 introduced by the 
differential rotation. 

According to the classification of 48] the above system of infinite coupled equations consists of two decoupled 
sub-systems among which there is no energy exchange. The first one is the so called "axial-led" system and consists 
of the polar radial I = perturbations that couple with the axial non-radial 1 = 1, the polar non-radial I = 2 and so 
on. The second system is the so called "polar-led" and consists of the polar dipole I — 1 perturbations that couple 
with the axial non-radial / = 2, the polar non-radial I — 3 and so on. Since the two sub-systems are decoupled in 
order to excite both we have to give independent appropriate initial data for time evolutions. In our study, for both 
polar and axial functions with I — 2 we assumed a generic gaussian pulse an an initial perturbation. The boundary 
conditions on the stellar surface are set by the vanishing of the Langrangian perturbation of the pressure AP = 0. 



IV. NUMERICAL RESULTS 



The perturbation equations derived in the previous sections have been studied numerically in order to compare the 
accuracy of the present approach against published results by nonlinear evolution codes [321 ]. More specifically we 
study the spectral properties of axisymmetric (m = 0) nonradial oscillations of a slowly and differentially rotating star, 
where the fluid is described by a barotropic EoS (T± = T). An extensive study of non-axisymmetric perturbations 
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of differentially rotating relativistic stars for a wide range of EoS will be the subject of a follow up paper [49j |. 
In the Cowling approximation, for barotropic and axisymmetric (m = 0) perturbations the system of evolution 
equations (|40| - ((44|) simplifies considerably and reduces to the following four coupled equations: 



Hi = 

uij = 

U2,l = 

U 3 ,l = 



-X' + 2v' \cl-v' 



e 2u 

ui.i + c s u 1J f e k ' - c s A—u 2 , 



2 f v 1 I (wi + 61373) — lu[ — 6lu' 3 



CI 



±1 



15 

T 



2 I V j TZ7 3 — CiJ 3 



Hi + j {2{w x + Gvu^Cf 1 - 15(fi 3 - 2u 3 )Cf 3 } u 3 j , 

-- {2(n7i + Gtua)^ 1 - 30m 3 Cf 3 } "2,i - ^ e~ 2A {((2 - 2rv') (w 1 + 6m 3 ) + r(zu[ + 6m' 3 )} C 

15 
T 



±1 

2 



[(2 - 2rv') w 3 + rvj' z ] Cf 3 \uij. 



(47) 
(48) 
(49) 

(50) 



where u 3 1 has been redefined as follows: 



u 3 ,i = u 3<t - -r 2 e 2l/ 



15 

(rui + 6w 3 )C 2 — —zu 3 £ 3 



Hi 



(51) 



Equation (|44[) of the variable £ r becomes trivial for barotropic pulsations, where Ti = V. In fact £ r is no longer an 
independent dynamical variable as one can argue from equations (|22[) and (|29p . 

The polar axisymmetric eigenfrequencies are symmetric with respect to the reversal of rotation direction, suggesting 
that the rotational corrections appear at second order in f2 [50j ]. Therefore, for first order slow rotation approximation 
the axisymmetric modes have the same frequencies as the nonrotating case. However, in equations l|47p - (|50p implicitly 
show up some 0(Vt 2 ) coupling terms, such as the last term of equations (|48|) -(|49 |) . In fact in accordance with 
equation (|50p . the axial velocity 1*3.; is a dynamical variable at 0(Q) perturbative order, while for a nonrotating star 
it is stationary and its profile can be chosen on the initial time slice. It is evident that this stationary initial condition 
cannot change the mode frequencies of polar perturbations. Here we use equations (|47 j) -(|50 |) without discarding the 
i*3 : 1 terms in equations (|4"8"| - (|49p . in order to investigate the 0(f2 2 ) corrections of the axial velocity on the polar mode 
frequencies. However, a complete second order slow rotation approximation is needed for a proper analysis that will 
take into account all coupling terms. 

The system of four coupled equations ([4T|) - ([5lI| has been evolved using second order accurate differencing methods 
both in time and in space based on the two step Lax-Wendroff algorithm. As a result, the numerical code achieved 
second order convergence. As initial conditions, for the enthalpy and the three velocity perturbations we set a Gaussian 
pulse, which excites several nonradial modes. The frequencies of these nonradial modes are then estimated by carrying 
out a Fast Fourier Transformation (FFT) of the resulting timeseries. 

In general, the system of perturbation equations has an infinite series of coupling terms, which depend on I. 
However, we got converging numerical results by keeping a finite number of equations, i.e. up to an I = Z max . In fact, 
starting from a simple configuration, where only a particular I is considered, one can gradually add up extra coupling 
equations until the results show convergence. An example is shown in table [III where we tabulated the frequencies of 
the fundamental quasi-radial mode (F) and the first overtone Hi for several simulations assuming different ^ max . The 
initial perturbed configuration is described by spherical oscillations, i.e. ^ max = 0. By gradually introducing the dipole 
l max = 1 and quadrupole ^ max = 2 oscillations, the quasi-radial mode frequencies have been improved. No further 
changes have been observed, when Z max = 3 (or higher) perturbations has been taken into account. Actually, in this 
case, convergence has been achieved already for Z max = 2. Since this behavior is shared also by the other nonradial 
modes, we restricted our study up to i max = 2. 

The accuracy of our description has been tested against published results from a nonlinear code (SAF from 
now on) for a sequence of stellar models. The frequencies of the first two quasi-radial (F and Hi) and two I = 2 
nonradial modes ( 2 / and 2 pi ) for some of the B-models are listed in table Hill Although these stellar models describe 
very fast rotating stars (see table [I]), the perturbative results show a good agreement (within the limits of the slow 
rotation approximation) with those obtained with nonlinear simulations in SAF. In fact, a difference of the order of 
7% appears only after the B6 model. Errors of this order are expected for very fast rotating stars since our results, 
which are based on the slow rotation approximation (first order in f2), cannot describe the flattening of the star due 
to the centrifugal force. 

An important aspect to emphasize is that our results do not show any splitting of the fundamental quasi-radial 
mode, as observed in SAF. This result supports the argument in [33|], that this effect (splitting) could be an artifact 
of the Cowling approximation in the non-linear regime. Moreover, we don't observe any significant change in the 
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Model 


F (kHz) 


Hi (kHz) 


V (kHz) 


2 Pi (kHz) 


BO 


2.706 (0%) 


4.547 (0%) 


1.846 (0%) 


4.100 (0%) 


Bl 


2.712 (2%) 


4.555 (2%) 


1.895 (<1%) 


4.117 (1%) 


B3 


2.735 (4%) 


4.578 (4%) 


1.915 (1%) 


4.124 (2%) 


B6 


2.797 (6%) 


4.624 (4%) 


1.944 (1%) 


4.134 (7%) 


B9 


2.885 (8%) 


4.686 (6%) 


1.974 (8%) 


4.147 (14%) 



TABLE III: Frequencies of the first two quasi-radial (F and Hi) and nonradial ( 2 f and 2 pi) modes, for differentially rotating 
B-models with A — 12 and Z max = 2. The relative difference between the perturbative and nonlinear results [32J is shown in 
parenthesis 




FIG. 5: The frequency of the 2 / mode as function of the parameter A is plotted for various B-models. For A = 500 km the star 
is practically rotates uniformly. The shifting of the mode frequency between the maximum and minimum values of A varies 
from 3% to 10% depending on the model. 



frequency of the 2 pi-mode in constrast to the results in SAF where an overall change of the order of 30% has been 
observed. 

The degree of differential rotation affects the mode frequencies due to the presence of extra coupling terms in the 
perturbation equations (|40|) - ((44|) . This effect can then be studied for several B-stellar models by varying the parameter 
A, from a highly differentially rotating configuration with A = 12 km to a uniformly rotating star with a very large 
A e.g. A = 500 km. As shown in figure IIVI the frequency of the nonradial fundamental 2 / mode increases as the 
parameter A increases for the whole sequence of B-models. In particular, the 2 / frequency already converges to its 
value for uniform rotation when A w 100 km. In this case, the angular velocity at the equator fl e is typically about 
2.5% smaller than fl c on the rotational axis. It is worth noticing that the effect of A on the mode frequencies is 
considerably larger for fast rotating stellar model. As shown in tabic HVl the same behavior has been observed also for 
the other quasi-radial and nonradial modes. Among these modes, the F-mode is the most sensitive to the differential 
rotation, as it can change up to 25% for the fastest rotating model B9. Instead, the 2 / and 2 pi modes can change up 
to 10 percent and 3 percent respectively. 
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Modes 


A [km] 


BO 


Bl 


B3 


B6 


B9 




12 


2.706 


2.702 


2.735 


2.797 


2.885 


F 


50 


2.706 


2.738 


2.849 


3.058 


3.352 




75 


2.706 


2.764 


2.954 


3.294 


3.747 




100 


2.706 


2.771 


2.960 


3.313 


3.777 




12 


1.846 


1.895 


1.915 


1.944 


1.974 


2 f 


50 


1.846 


1.928 


1.987 


2.058 


2.117 




75 


1.846 


1.947 


2.038 


2.111 


2.169 




100 


1.846 


1.954 


2.040 


2.118 


2.170 




12 


4.100 


4.117 


4.124 


4.134 


4.147 




50 


4.100 


4.124 


4.143 


4.163 


4.196 




75 


4.100 


4.130 


4.155 


4.175 


4.209 




100 


4.100 


4.134 


4.157 


4.176 


4.210 



TABLE IV: Frequencies, in kHz, of the quasi-radial F mode, and 2 / and 2 pi nonradial modes for B-stellar models with different 
rotational parameter A. 



V. CONCLUSIONS AND DISCUSSION 

In this article we derived the general equations describing the perturbations of slowly and differentially rotating 
neutron stars in the Cowling approximation. The equations have been derived for spherical stars (slow rotation) on 
which the differential rotation has been described via the j-constant law. 

As a first step and test, results have been derived for axisymmetric, barotropic oscillations. In the system of 
perturbation equations some terms of second order 0(fl 2 ) have been maintained and their effects on the axisymmetric 
eigenfrequencies has been studied. The comparison of the derived results with those produced by the full non-linear 
code [32;] suggests that this perturbative approach can provide quite accurate estimates of the oscillation frequencies 
within the limits of the slow rotation approximation. A point to be stressed is that since we keep only terms first 
order in rotation 0(fl) the background stellar model is still spherical, and the compactness of the star M/R remains 
the same along the sequence of stellar models. This seems to be one of the most important reasons for the deviations 
observed between our results and those of SAF for the faster rotating members of the B-model sequence. It is expected 
that the inclusion of second order terms in rotation, which will take into account the flattening of the star, will improve 
the agreement between the perturbative results with the non-linear ones for the very fast rotating models. 

An advantage of this perturbative approach is the numerical simplicity and the ability for an easier physical 
interpretation of the derived results. In a follow up paper we study the non-axisymmetric oscillations of differentially 
rotating relativistic stars for various EoS with special emphasis to the hot EoS. 
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APPENDIX A: HARMONIC COMPONENTS OF THE J-CONSTANT ROTATION LAW 

In this Appendix we present the explicit expression for the terms Xf and X\ described in equations (fl"T)|) and (|13[) . 
Since the series expansions ©-(JT]) for to and Q have been truncated for I > 3, and that terms with even I vanish due 
to the rotation law parity, we show the expressions for harmonic indices / = 1 and 1 = 3. The two nearly "Newtonian 
terms" Xf are: 



N _ 

1 ~ 2 x 2 



N _ ]_A_ 

±3 ~ 12 x 2 C 



A 2 



In 



yj ^A 2 + x 2 + X 



X^A 2 + x 2 \^A 2 + x 2 
.A 2 6 A 2 



X, 



sjA 2 + X 2 +X 



5 A 2 \ 

1 + 7.5'— (1 + -— In, 

x 2 xV A2 + x 2 V 4x 2 / \ J^ A 2 + x 2~ / 



The four terms that contribute in the estimation of the "relativistic corrections" Xr are: 



Ji 



2 X 2 - 3A 2 



: m 1 



V^ 2 + x 2 V^v^TF-x, 



.4 2 



3 7T 

7^ 
x In 



x^ 2 + x 2 VvV^+r^-x 

A 4 



-525yl b - 490AV - 56AV + 16x 



+ X 2 



(525A 4 + 840A 2 x 2 + 336x 4 



V^ 2 + X 2 + X 

,x/^4 2 + x 2 -x, 

where we have set x = ^e - ". 



APPENDIX B: COEFFICIENTS OF THE PERTURBED CONSERVATION EQUATIONS 

In this Appendix we present the coefficients of the perturbation equations (|3"rj)) -(|3"5 |) . 
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7 W 

'4,1m 



ai.im = -U2,i m +Hi m -im(Cli+6fl 3 )u2,im, (B9) 

02, im = 2 (wi + 6t37 3 )u 3i ; m , (BIO) 

15 

a3,im = y imft 3 u 2j im , (Bll) 

4jm = -15(fi3-2w 3 ) u 3 , im , (B12) 

-30tu 3 u 3iim , (B13) 

h,im = -v>3,im -im(f2i +6fi 3 )u3,im (B14) 

&2,im = -2 [tUi + 6tU 3 ] M 2 , im , (B15) 

15 

&3,im = y«™fi 3 U 3j ; m , (B16) 
& Sm = 15(^3-2^3) U2,lm, (B17) 

& i? m = 30n7 3M2 , im , (B18) 

c Um = c l ( ro i + 6 ^3) (r 2 e~ 2A lm - Au 2 , lm) 

- re~ 2X { [2 - rv' - c 2 (2 + r (2i/ - A'))] {vJt + 6 m 3 ) + r [w[ + 6 zu' 3 )} Ui, im 

= r 2 e~ 2v {w 1 + 6w 3 ) H lm - re' 2 * [(2 - 2rv'){w x + 6w 3 ) + r{w[ + 6*4)] Ui,; m , (B19) 

15 g / 1 9 — 2A / \ 

c 2,/m = yC s w 3 (A«2,im-re u Mm ) 

+ yre- 2A {[2-r^- C 2(2 + r(2^-A'))] 073 + rzj' 3 } u 1Jm 

= -^-r 2 e- 2v w 3 Hi m + ^-re- 2X [{2-2rv')w 3 +rm' 3 ]u 1 ,i m , (B20) 



where we have used the definitions, 

tui(r) = Qi(r) - wi(r) , 
tu 3 (r) = fi 3 (r) - w 3 (r) . 

APPENDIX C: ANGULAR INTEGRATION 

Here, we define the linear operators C i J with i, j S N, which come out from the angular integration of equations (|36|) - 
(|39[) . This procedure relies on the orthogonality properties of the spherical harmonics, on the following angular 
integrals: 



(dVp m ,dY lm , 1 OY^dY, 



52 \ 96* 96* sin 90 90 



S'2 



*l^ + ^ii.)sf-fc«-, «s> 



and the two relations: 



COS0Y; m — Ql + lmYl + lm + QlmXl-\m. , (C3) 

sin6»9 e y /m = Qi+i m /F i+ i m - Q im (Z + 1) Y/_i m , (C4) 
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where Qi m is defined as follows: 



/ (Z - m) (Z + m) 
Qlm ~\j {2l-l){2l + l)- 

The operators £^ bl which introduce couplings of a perturbation term/expression with harmonic index Z with the 
term/expression having harmonic indices Z ± 1 are given by the following relations: 



(C6) 
(C7) 



C^A = V A Vm < \ dCl Y{ m sin d Y Vm , = (I - 1) Qi m Ai- lm - (I + 2) Qi+i m A l+lm , 

v, JS 2 
I'm' 

d^A = V Ai'm' / d& d e Y* m sin 6 Y Vm > = - (I + 1) Qi m Ai-im + lQi+i m Ai+i m , 

I'm' J $ 2 

£±U = y^Ai'm' I dCl (dgY^dgYw + -L-d^d^YvJ] COS0 
ft' J s* V sin 6 J 

= G — 1) G + 1) Qlm-A-l-lm +1(1 + 2) Ql + im-A-i + im , (C8) 
£4* .A EE V^l/'m' / dfi COS Yi/ m / = Qim-Ai-lm + Q(+lm»4(+lm , (C9) 

£±U = V A' m ' / dn (^y^ ^y,/ m / - fyy^ fl y,, m , ) sine - im (rf 1 + £± x ) . 



(CIO) 



Notice that the operator £4 can be defined in terms of the £f and Cf operators as follows: 

£ 4 ±1 = -^(£r 1 + A ±1 ) 

In a similar way the operators £* 2 that introduce couplings of a perturbation expression with harmonic index I with 
term/expression with harmonic indices / ± 2 are: 

Cf 2 A = V A Vm ' \ dflY^ sm 2 9Y llm , = -Qi-i m Qi m Ai- 2 m + (l - <?L - Q?+lm) 

-QlmQl + lmAl+2m , (Gil) 

£± 2 „4 EE Y^'m' / dfi^yim Bind COS 0y Z , m - = - (/ + 1) Qi-l TO Ql m ^J_2m + [lQ 2 +lm - (I + 1) Qf m ] Aim 

+'Qi+lm<3i+2m-4i+2m , (C12) 

Cf 2 A ee Y^Avm' f dCl Y\ m sin 6 cos 0$Yi' m / = (I - 2) Q,_i m Q Im ^_ 2m + [*Q 2 +i m - (/ + 1)QL] -4™ 

I'm' J S 2 

— (I + 3) Ql+lmQl+2mAl+2m , (C13) 

£± 2 .a ee V ^, m , / rfo f ^y*, a e y, m , + d^Y Vm \ sin 2 = - (i - 2) (i + 1) Qi-i m Qi m Ai- 2m 

w V sin 6 J 

+ [l(l + l)-(l + 1) (I - 2) Q 2 m - I (I + 3) Q 2 +lm ] Aim -1(1 + 3) Ql + lmQl+2mAl+2m ■ (C14) 

Finally, the operators £^ 3 that introduce couplings of a perturbation expression with harmonic index Z with 
terms/expressions with harmonic indices I ± 3 are: 
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Cf 3 A = y^A'm' I dClY^ sin 3 doY Vm , 
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Cf 3 A 



+Qhn [iQl-lm + _ 1) (l _ Qlm ~ Ql+lm)] A-lm 
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(C15) 
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I'm' 6 



/ — 3m 



1-lm 
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cZf2 <9 e Y^ sin 3 F z - m / = (Z + 1) Q\-imQi-\ m QimA 



I — 3m 



-Q im [(J + l) + ZQ 2 +lm - (I + 1) (Q?_ lm + QL)] A-u 

+Ql+lm [l + {l + 1) QL - I (Qf +lm + Qf+2 J] ^l+lm 
— ZQi + lmQz+2mQi+3m-4z+3m • 



(C16) 



(C17) 



APPENDIX D: ROTATION LAW IN ISOTROPIC COORDINATES 

In this Appendix, we describe the j-constant rotation law in isotropic coordinates, and the appropriate transforma- 
tion in Schwarzshild coordinates. In isotropic coordinates the axially symmetric spacetime for a slowly rotating body 
can be described by the following metric: 

ds 2 = -e 2a dt 2 + e 2p [df 2 + f 2 (d9 2 + sin 2 6d(f> 2 )] -2e 2p uj f 2 sin 2 dt d<P , (Dl) 



where a and (3 are functions of the isotropic coordinate f. Equation (|Dip has been derived by the slow rotating 
metric ((T|) in Schwarschild coordinate via the following relations: 

dr r \ r ) 

e = t, (D3) 
r 

where M = M(r) is the mass function in Schwarzshild coordinates. At the stellar surface the metric functions must 
smoothly join the Schwarzshild solution, this leads to the following matching conditions in isotropic coordinates 

1 - — ) ( 1 + — ) > (° 4 ) 

J> - 



2f j \ 2f 

M N 2 



1 + Yf) • ^ 

The stellar radius R 7 in isotropic coordinates, is given in terms of the stellar mass M and the radius R in Schwarzshild 
coordinates by the following relation: 

R=-(r-M+ \/R 2 - 2MR) . (D6) 

Finally, the relativistic j-constant rotation law can be determined via the relation: 

= A 2 {Q c - Q) , (D7) 
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which leads to the following expression of the angular velocity profile: 

A 2 tt c + e 2 ^-") f 2 sin 2 9u(f,0) 



m0) = V (D8) 



^2 _|_ e 2((3-a) p 



sm 



where A is a parameter that describes the degree of differential rotation in isotropic coordinates. By defining a new 
function 

x = e^ a 4 sin 6», (D9) 
A 

we can rewrite equation (|D8|) in a more compact form i.e. 

fi(f,e) _ 1 + x 2 uj(f,e)/n c 



tt c 1 + X 2 



(D10) 



It is worth noticing that equation ([5]) can be written in the same form of equation (|D10[) if the function x is replaced 
by 

x = e~ v — sin 9 . (Dll) 

Equations (|D8|) and ([5]) are obviously related by the coordinate transformation between isotropic and Schwarzshild 
coordinates. A typical choice for A is to set its value equal to the stellar radius R, which, as we show later, corresponds 
to f2 e ~ f2 c /3. This relation between the angular velocity at the equator and the rotation axis seems to be in good 
agreement with the rotation patters estimated in the remnants of hypermassive binary mergers [J]. The matching 
conditions (|D4[) and (|D5[) lead to the following expression for x e at the equator: 

M\ 3 A M \ 1 R 



, \ ■ 1 (D12) 

2RJ V 2RJ A 

A coordinate transformation leads to the following expression for x e in Schwarzshild coordinates: 

2M\ _1/2 R 



x*=(l--n) J- (D13) 

From equations (|D10|) - (|D13|) . it is worth noticing that the ratio VL e /Sl c depends practically on the compactness of 
the star. In order to set up an equilibrium configuration in Schwarzshild coordinates, which has the same value of 
il e /il c as the one used in non-linear calculations in isotropic coordinates, the natural choice is A = A as it comes out 
from the equations (|D13p and (|D11[) . For example, for a polytropic stellar model with mass M = 1.4 M© and radius 
R = 14.151 km, equation (|D6[) gives R = 12 km and leads to i e = 1.4. By neglecting the corrections given by the 
metric function u> in equations (|D8|) and ([5|), the choice A = R leads to x e = x e and tt e = 0.33 fi c in both coordinate 
systems. 
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